Project

Mapping of Ki67 interactions with the genome and comparison with lamina interactions.

Introduction

Various analyses of HCT116 cells before and after Ki67 loss using a Ki67-AID degron system.

Method

NA

Set-up

Set the parameters and list the data.

# Load dependencies - without warnings / messages
library(tidyverse)
library(GenomicRanges)
library(rtracklayer)
library(ggplot2)
library(RColorBrewer)
library(GGally)
library(corrr)
library(caTools)
library(ggrastr)
library(colorspace)

# Prepare output 
output_dir <- "ts220503_7_hct116_ki67_aid"
dir.create(output_dir, showWarnings = FALSE)


# Load input
chromosomes <- c(paste0("chr", 1:22), "chrX")


input_dir <- "ts220503_1_data_gathering"

bin_size <- readRDS(file.path(input_dir, "bin_size.rds"))
centromeres <- readRDS(file.path(input_dir, "centromeres.rds"))

colors_set1 <- readRDS(file.path(input_dir, "colors_set1.rds"))
colors_set2 <- readRDS(file.path(input_dir, "colors_set2.rds"))
colors_set3 <- readRDS(file.path(input_dir, "colors_set3.rds"))

tib_padamid_replicates <- readRDS(file.path(input_dir, 
                                            "tib_padamid_replicates.rds"))
tib_padamid_combined <- readRDS(file.path(input_dir, 
                                          "tib_padamid_combined.rds"))

gr_padamid_replicates <- readRDS(file.path(input_dir, 
                                           "gr_padamid_replicates.rds"))
gr_padamid_combined <- readRDS(file.path(input_dir, 
                                         "gr_padamid_combined.rds"))

tib_hmm_replicates <- readRDS(file.path(input_dir, "tib_hmm_replicates.rds"))
tib_hmm_combined <- readRDS(file.path(input_dir, "tib_hmm_combined.rds"))

padamid_metadata_replicates <- readRDS(file.path(input_dir, 
                                                 "padamid_metadata_replicates.rds"))
padamid_metadata_combined <- readRDS(file.path(input_dir, 
                                               "padamid_metadata_combined.rds"))



# Prepare seqnames
chrom_sizes <- tibble(seqnames = seqlevels(gr_padamid_combined),
                      length = seqlengths(gr_padamid_combined)) %>%
  arrange(-length)


# Scale pA-DamID scores?
tib_padamid_combined <- tib_padamid_combined %>%
  mutate_at(4:ncol(.), function(x) scale(x)[, 1]) %>%
  filter(seqnames != "chrY")
library(knitr)
opts_chunk$set(fig.width = 10, fig.height = 4, cache = T,
               dev=c('png', 'pdf'), fig.path = file.path(output_dir, "figures/")) 
pdf.options(useDingbats = FALSE)
# From Fede:
# ggpairs custom functions
corColor <- function(data, mapping, color = I("black"), sizeRange = c(1, 3), ...) {

  x   <- eval_data_col(data, mapping$x)
  y   <- eval_data_col(data, mapping$y)
  r   <- cor(x, y)
  rt  <- format(r, digits = 3)
  tt  <- as.character(rt)
  cex <- max(sizeRange)

  # helper function to calculate a useable size
  percent_of_range <- function(percent, range) {
    percent * diff(range) + min(range, na.rm = TRUE)
  }

  # plot correlation coefficient
  p <- ggally_text(label = tt, mapping = aes(), xP = 0.5, yP = 0.5,
                   # size = I(percent_of_range(cex * abs(r), sizeRange)), 
                   size = 6, 
                   color = color, ...) +
    theme(panel.grid.minor=element_blank(),
          panel.grid.major=element_blank())

  corColors <- RColorBrewer::brewer.pal(n = 7, name = "RdYlBu")[2:6]

  if (r <= boundaries[1]) {
    corCol <- corColors[1]
  } else if (r <= boundaries[2]) {
    corCol <- corColors[2]
  } else if (r < boundaries[3]) {
    corCol <- corColors[3]
  } else if (r < boundaries[4]) {
    corCol <- corColors[4]
  } else {
    corCol <- corColors[5]
  }

  p <- p +
    theme(panel.background = element_rect(fill = corCol))

  return(p)
}

customScatter <- function (data, mapping) 
{
    p <- ggplot(data = data, mapping) + 
      geom_bin2d(bins = 100) +
      geom_smooth(method = "lm", se = T, col = "red") +
      scale_fill_gradient(low = "lightgrey", high = "black", name = "Count") +
      theme_bw()
    
    p 
}

PlotScatter <- function(tib, n1, n2, color_by = NULL, identity = F, smooth = 1) {
  # Get tibble
  tib <- tib %>%
    dplyr::select(seqnames, matches(n1), matches(n2)) %>%
    rename_all(~ c("seqnames", "n1", "n2"))
  
  if (smooth > 1) {
    tib <- tib %>%
      group_by(seqnames) %>%
      mutate(n1 = runmean(n1, k = smooth),
             n2 = runmean(n2, k = smooth))
  }
  
  # Prepare color
  if (! is.null(color_by)) {
    tib <- tib %>%
      add_column(color = color_by) %>%
      drop_na()
    alpha = 1
    limits_color <- quantile(tib$color, c(0.001, 0.999), na.rm = T)
    tib$color[tib$color < limits_color[1]] <- limits_color[1]
    tib$color[tib$color > limits_color[2]] <- limits_color[2]
  } else {
    tib <- tib %>% drop_na()
    tib$color = "1"
    alpha = 0.02
  }
  
  # Plot
  xlimits <- quantile(tib$n1, c(0.001, 0.999), na.rm = T) * 1.4
  ylimits <- quantile(tib$n2, c(0.001, 0.999), na.rm = T) * 1.4
  
  plt <- tib %>%
    arrange(sample(1:nrow(.), size = nrow(.), replace = F)) %>%
    ggplot(aes(x = n1, y = n2, color = color)) +
    geom_point(size = 0.5, alpha = alpha) +
    geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
    geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
    xlab(n1) +
    ylab(n2) +
    ggtitle(paste0("Spearman: ", 
                   round(cor(tib$n1, tib$n2, use = "complete.obs",
                             method = "spearman"), 2))) +
    coord_cartesian(xlim = xlimits, ylim = ylimits) +
    theme_bw() +
    theme(aspect.ratio = 1)
  
  # Prepare color
  if (! is.null(color_by)) {
    plt <- plt +
      scale_color_gradient2(low = "blue", mid = "grey", high = "red",
                            midpoint = 0)
  } else {
    plt <- plt + 
      scale_color_manual(values = "black", guide = F)
  }
  if (identity) plt <- plt + geom_abline(slope = 1, intercept = 0, col = "red", linetype = "dashed")
  
  plot(plt)
  
}

PlotDataTracksLight <- function(input_tib, exp_names, centromeres, 
                                color_groups, plot_chr = "chr1", 
                                plot_start = 1, plot_end = 500e6,
                                smooth = 1, color_list = NULL,
                                fix_yaxis = F, aspect_ratio = NULL,
                                lighten_negative = T, raster = T,
                                ylimits = NULL) {
  
  # Get the scores for the samples
  tib_plot <- input_tib %>%
    dplyr::select(seqnames, start, end, 
                  all_of(exp_names))
  
  if (smooth > 1) {
    tib_plot <- tib_plot %>%
      mutate_at(vars(contains("_")), 
                runmean, k = smooth)
  }
  
  # Filter for plotting window
  tib_plot <- tib_plot %>%
    filter(seqnames == plot_chr,
           start >= plot_start,
           end <= plot_end)
  
  # Gather
  tib_gather <- tib_plot %>%
    gather(key, value, -seqnames, -start, -end) %>%
    mutate(key = factor(key, levels = exp_names),
           fill_column = color_groups[match(key, names(color_groups))],
           fill_column = factor(fill_column, levels = unique(color_groups)))
  
  # If desired, make negative values a lighter shade to improve visualization
  if (lighten_negative) {
    tib_gather <- tib_gather %>%
      mutate(fill_column = interaction(fill_column,
                                       value < 0))
  }
  
  
  # Centromeres
  centromeres.plt <- as_tibble(centromeres) %>%
    filter(seqnames == plot_chr) %>%
    gather(key_centromeres, value, start, end)
  
  # Plot
  if (is.null(ylimits)) {
    ylimits <- quantile(tib_gather$value, c(0.001, 0.999), na.rm = T)
  }
  
  fill_column <- NULL
  
  plt <- tib_gather %>% 
    ggplot(aes(fill = fill_column))
  
  
  # Set all counts tracks to the same limits
  if (raster) {
    plt <- plt + 
      rasterize(geom_rect(aes(xmin = start / 1e6, xmax = end / 1e6, 
                              ymin = 0, ymax = value)),
                dpi = 300)
  } else {
    plt <- plt + 
      geom_rect(aes(xmin = start / 1e6, xmax = end / 1e6, 
                    ymin = 0, ymax = value))
  }
  
  plt <- plt + 
    geom_hline(yintercept = 0, size = 0.5) +
    geom_line(data = centromeres.plt, 
              aes(x = value / 1e6, y = 0),
              color = "black", size = 2) +
    facet_grid(key ~ ., scales = "free_y") +
    xlab(paste0(plot_chr, " (Mb)")) +
    ylab("Score") +
    scale_x_continuous(expand = c(0, 0)) + 
    scale_y_continuous(expand = c(0.025, 0.025)) +
    theme_classic()
  
  if (! is.null(color_list)) {
    
    if (lighten_negative) {
      color_list <- c(color_list,
                      lighten(color_list, amount = 0.5))
    }
    
    colors <- levels(tib_gather$fill_column)
    
    color_list <- color_list[1:length(colors)]
    names(color_list) <- colors
    #colors <- recode(colors, !!!color_list)
    
    plt <- plt +
      scale_fill_manual(values = color_list, guide = "none")
  } else {
    plt <- plt +
      scale_fill_brewer(palette = "Set1", guide = "none")
  }
  
  if (fix_yaxis) {
    plt <- plt + 
      coord_cartesian(xlim = c(max(c(plot_start,
                                     min(tib_gather$start))) / 1e6,
                               min(c(plot_end,
                                     max(tib_gather$end))) / 1e6),
                      ylim = ylimits)
  } else {
    plt <- plt + 
      coord_cartesian(xlim = c(max(c(plot_start,
                                     min(tib_gather$start))) / 1e6,
                               min(c(plot_end,
                                     max(tib_gather$end))) / 1e6))
  }
  
  if (! is.null(aspect_ratio)) {
    plt <- plt +
      theme(aspect.ratio = aspect_ratio)
  }
  
  plot(plt)
  
}

# Similar data tracks as the function above, but with the functionality to plot
# matching data tracks on top of each other
PlotDataTracksLightFaceted <- function(input_tib, exp_names, centromeres, 
                                       color_groups, facet_groups,
                                       plot_chr = "chr1", size = 0.75,
                                       plot_start = 1, plot_end = 500e6,
                                       smooth = 1, color_list = NULL,
                                       fix_yaxis = F, scale_tracks = F, 
                                       aspect_ratio = NULL, raster = T, ylimits = NULL) {
  
  # # Exp names is a vector of sample names
  # exp_search <- paste(exp_names, collapse = "|")
  
  # Get the scores for the samples
  tib_plot <- input_tib %>%
    dplyr::select(seqnames, start, end, 
                  all_of(exp_names))
  
  if (smooth > 1) {
    tib_plot <- tib_plot %>%
      mutate_at(vars(contains("_")), 
                runmean, k = smooth)
  }
  
  if (scale_tracks) {
    tib_plot <- tib_plot %>%
      mutate_at(vars(contains("_")), 
                function(x) scale(x)[, 1])
  }
  
  # Filter for plotting window
  tib_plot <- tib_plot %>%
    filter(seqnames == plot_chr,
           start >= plot_start,
           end <= plot_end)
  
  # Gather
  tib_gather <- tib_plot %>%
    gather(key, value, -seqnames, -start, -end) %>%
    mutate(key = factor(key, levels = exp_names),
           fill_column = color_groups[match(key, names(color_groups))],
           fill_column = factor(fill_column, levels = unique(color_groups)),
           facet_column = facet_groups[match(key, names(facet_groups))],
           facet_column = factor(facet_column, levels = unique(facet_groups)))
  
  
  # Centromeres
  centromeres.plt <- as_tibble(centromeres) %>%
    filter(seqnames == plot_chr) %>%
    gather(key_centromeres, value, start, end)
  
  # Plot
  if (is.null(ylimits)) {
    ylimits <- quantile(tib_gather$value, c(0.001, 0.999), na.rm = T)
  }
  
  fill_column <- NULL
  
  plt <- tib_gather %>% 
    ggplot(aes(fill = fill_column))
  
  
  # Set all counts tracks to the same limits
  plt <- plt + 
    #geom_rect(aes(xmin = start / 1e6, xmax = end / 1e6, 
    #              ymin = 0, ymax = value)) +
    geom_hline(yintercept = 0, size = 0.5) +
    geom_line(data = centromeres.plt, 
              aes(x = value / 1e6, y = 0),
              color = "black", size = 2) +
    facet_grid(facet_column ~ ., scales = "free_y") +
    xlab(paste0(plot_chr, " (Mb)")) +
    ylab("Score") +
    scale_x_continuous(expand = c(0, 0)) + 
    scale_y_continuous(expand = c(0.025, 0.025)) +
    theme_classic()
  
  if (raster) {
    plt <- plt +
      rasterize(geom_line(aes(x = (start+end)/2e6, y = value,
                              col = fill_column), size = size), 
                dpi = 300)
  } else {
    plt <- plt +
      geom_line(aes(x = (start+end)/2e6, y = value,
                  col = fill_column), size = size)
  }
  
  if (! is.null(color_list)) {
    colors <- levels(tib_gather$fill_column)
    
    color_list <- color_list[1:length(colors)]
    names(color_list) <- colors
    #colors <- recode(colors, !!!color_list)
    
    plt <- plt +
      scale_fill_manual(values = color_list, guide = "none") +
      scale_color_manual(values = color_list, guide = "none")
  } else {
    plt <- plt +
      scale_fill_brewer(palette = "Set1", guide = "none") +
      scale_color_brewer(palette = "Set1", guide = "none")
  }
  
  if (fix_yaxis) {
    plt <- plt + 
      coord_cartesian(xlim = c(max(c(plot_start,
                                     min(tib_gather$start))) / 1e6,
                               min(c(plot_end,
                                     max(tib_gather$end))) / 1e6),
                      ylim = ylimits)
  } else {
    plt <- plt + 
      coord_cartesian(xlim = c(max(c(plot_start,
                                     min(tib_gather$start))) / 1e6,
                               min(c(plot_end,
                                     max(tib_gather$end))) / 1e6))
  }
  
  if (! is.null(aspect_ratio)) {
    plt <- plt +
      theme(aspect.ratio = aspect_ratio)
  }
  
  plot(plt)
  
}

PlotScatterBinned <- function(tib, n1, n2, color_by = NULL, identity = F, 
                              n_min = 10, ylimits_col = c(-2.4, 2.4),
                              count_range = c(0, 400), smooth = 1,
                              remove_outliers = F) {
  # Get tibble
  tib_process <- tib %>%
    dplyr::select(seqnames, all_of(n1), all_of(n2)) %>%
    rename_all(~ c("seqnames", "n1", "n2"))
  
  # Calculate pearson correlation without any smoothing
  tib_cor <- round(cor(tib_process$n1, tib_process$n2, use = "complete.obs",
                       method = "pearson"), 2)
  
  if (smooth > 1) {
    tib_process <- tib_process %>%
      group_by(seqnames) %>%
      mutate(n1 = runmean(n1, k = smooth),
             n2 = runmean(n2, k = smooth))
  }
  
  
  if (! is.null(color_by)) {
    tib_process <- tib_process %>%
      add_column(color = color_by)
  }
  
  tib_process <- tib_process %>%
    drop_na()
  
  # Change color range
  if (! is.null(color_by)) {
    limits_color <- quantile(tib_process$color, c(0.001, 0.999), na.rm = T)
    tib_process$color[tib_process$color < limits_color[1]] <- limits_color[1]
    tib_process$color[tib_process$color > limits_color[2]] <- limits_color[2]
  }
  
  # Replace outliers
  if (remove_outliers) {
    tib_process <- tib_process %>%
      mutate(n1 = case_when(n1 > quantile(n1, 0.999) ~ quantile(n1, 0.999),
                            n1 < quantile(n1, 0.001) ~ quantile(n1, 0.001),
                            T ~ n1),
             n2 = case_when(n2 > quantile(n2, 0.999) ~ quantile(n2, 0.999),
                            n2 < quantile(n2, 0.001) ~ quantile(n2, 0.001),
                            T ~ n2))
  }
  
  # Metrics
  n1_min = min(tib_process$n1) - 0.001
  n1_max = max(tib_process$n1) + 0.001
  n1_binsize <- (n1_max - n1_min) / 49
  
  n2_min = min(tib_process$n2) - 0.001
  n2_max = max(tib_process$n2) + 0.001
  n2_binsize <- (n2_max - n2_min) / 49
  
  tib_summary <- tib_process %>%
    mutate(n1_cut = cut(n1, 
                        seq(n1_min, n1_max, length.out = 50)),
           n2_cut = cut(n2, 
                        seq(n2_min, n2_max, length.out = 50))) %>%
    mutate(n1_bin = as.numeric(as.factor(n1_cut)),
           n2_bin = as.numeric(as.factor(n2_cut))) %>%
    mutate(n1_bin = n1_min - n1_binsize/2 + n1_bin * n1_binsize,
           n2_bin = n2_min - n2_binsize/2 + n2_bin * n2_binsize) %>%
    group_by(n1_bin, n2_bin)
  
  # Plot
  if (! is.null(color_by)) {
    tib_summary <- tib_summary %>%
    dplyr::summarise(n = n(),
                     mark = mean(color)) %>%
    ungroup() %>%
    filter(n >= n_min)
    
    plt <- tib_summary %>%
      ggplot(aes(x = n1_bin, y = n2_bin)) +
      geom_tile(aes(fill = mark)) +
      scale_fill_gradient2(low = "blue", mid = "grey", high = "red",
                           midpoint = 0, limits = ylimits_col, 
                           na.value = "green")
  } else {
    tib_summary <- tib_summary %>%
    dplyr::summarise(n = n()) %>%
    ungroup() %>%
    filter(n >= n_min)
    
    plt <- tib_summary %>%
      ggplot(aes(x = n1_bin, y = n2_bin)) +
      geom_tile(aes(fill = n)) +
      scale_fill_gradient(low = "lightgrey", high = "black", name = "Count",
                          limits = count_range, na.value = "green")
  }
  
  plt <- plt + 
    geom_hline(yintercept = 0, linetype = "dashed", col = "black") +
    geom_vline(xintercept = 0, linetype = "dashed", col = "black") +
    xlab(n1) +
    ylab(n2) +
    ggtitle(paste0("Pearson: ", tib_cor)) +
    theme_bw() +
    theme(aspect.ratio = 1)
  
  if (identity) plt <- plt + geom_abline(slope = 1, intercept = 0, 
                                         col = "black", linetype = "dashed")
  
  plot(plt)
  
}

1. Scatter plots for Ki67 depletion

Let’s start with making some scatter plots.

# Get metadata
padamid_metadata_combined_aid <- padamid_metadata_combined %>%
  filter(experiment == "Ki67_aid") %>%
  filter(! grepl("control", condition),
         ! grepl("thymidine", condition))

# Get tibble
tib <- tib_padamid_combined %>%
  mutate(# RO vs double-thy
         RO_vs_doublethy_Ki67 = HCT116_AID_RO_Ki67 - HCT116_AID_doublethy_Ki67,
         RO_vs_doublethy_LMNB1 = HCT116_AID_RO_LMNB1 - HCT116_AID_doublethy_LMNB1,
         RO_vs_doublethy_H3K27me3 = HCT116_AID_RO_H3K27me3 - HCT116_AID_doublethy_H3K27me3,
         RO_vs_doublethy_H3K9me3 = HCT116_AID_RO_H3K9me3 - HCT116_AID_doublethy_H3K9me3,
         # IAA vs nonIAA, double-thy
         doublethy_LMNB1_diff = HCT116_AID_doublethy_IAA_LMNB1 - HCT116_AID_doublethy_LMNB1,
         doublethy_H3K27me3_diff = HCT116_AID_doublethy_IAA_H3K27me3 - HCT116_AID_doublethy_H3K27me3,
         doublethy_H3K9me3_diff = HCT116_AID_doublethy_IAA_H3K9me3 - HCT116_AID_doublethy_H3K9me3,
         # IAA vs nonIAA, RO
         RO_LMNB1_diff = HCT116_AID_RO_IAA_LMNB1 - HCT116_AID_RO_LMNB1,
         RO_H3K27me3_diff = HCT116_AID_RO_IAA_H3K27me3 - HCT116_AID_RO_H3K27me3,
         RO_H3K9me3_diff = HCT116_AID_RO_IAA_H3K9me3 - HCT116_AID_RO_H3K9me3)

# Smoothing of differences?
smooth_differences = F

if (smooth_differences) {
  tib <- tib %>%
    group_by(seqnames) %>%
    mutate_at(vars(contains("vs"), contains("diff")), 
              function(x) runmean(x, k = 3)) %>%
    ungroup()
}




# Simple correlations
# 1) RO vs doublethy
PlotScatterBinned(tib, "HCT116_AID_doublethy_Ki67", "HCT116_AID_RO_Ki67", 
                  identity = T, count_range = c(0, 800))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "HCT116_AID_doublethy_LMNB1", "HCT116_AID_RO_LMNB1", 
                  identity = T, count_range = c(0, 800))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "HCT116_AID_doublethy_H3K27me3", "HCT116_AID_RO_H3K27me3", 
                  identity = T, count_range = c(0, 800))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "HCT116_AID_doublethy_H3K9me3", "HCT116_AID_RO_H3K9me3", 
                  identity = T, count_range = c(0, 800))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# 2) Effect of IAA
PlotScatterBinned(tib, "HCT116_AID_doublethy_LMNB1", "HCT116_AID_doublethy_IAA_LMNB1", 
                  identity = T, count_range = c(0, 800))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "HCT116_AID_doublethy_H3K27me3", "HCT116_AID_doublethy_IAA_H3K27me3", 
                  identity = T, count_range = c(0, 800))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "HCT116_AID_doublethy_H3K9me3", "HCT116_AID_doublethy_IAA_H3K9me3", 
                  identity = T, count_range = c(0, 800))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "HCT116_AID_RO_LMNB1", "HCT116_AID_RO_IAA_LMNB1", 
                  identity = T, count_range = c(0, 800))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "HCT116_AID_RO_H3K27me3", "HCT116_AID_RO_IAA_H3K27me3", 
                  identity = T, count_range = c(0, 800))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, "HCT116_AID_RO_H3K9me3", "HCT116_AID_RO_IAA_H3K9me3", 
                  identity = T, count_range = c(0, 800))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# Simple data tracks - for double-thy & RO
PlotDataTracksLight(input_tib = tib, 
                    exp_names = c("HCT116_AID_doublethy_Ki67",
                                  "HCT116_AID_doublethy_LMNB1",
                                  "HCT116_AID_doublethy_IAA_LMNB1",
                                  "HCT116_AID_doublethy_H3K27me3",
                                  "HCT116_AID_doublethy_IAA_H3K27me3",
                                  "HCT116_AID_doublethy_H3K9me3",
                                  "HCT116_AID_doublethy_IAA_H3K9me3"),
                    centromeres,
                    color_groups = c(HCT116_AID_doublethy_Ki67 = "Ki67",
                                     HCT116_AID_doublethy_LMNB1 = "LMNB1",
                                     HCT116_AID_doublethy_IAA_LMNB1 = "LMNB1",
                                     HCT116_AID_doublethy_H3K27me3 = "H3K27me3",
                                     HCT116_AID_doublethy_IAA_H3K27me3 = "H3K27me3",
                                     HCT116_AID_doublethy_H3K9me3 = "H3K9me3",
                                     HCT116_AID_doublethy_IAA_H3K9me3 = "H3K9me3"),
                    smooth = 9, plot_chr = "chr11", fix_yaxis = T,
                    plot_start = 25e6, plot_end = 100e6,
                    color_list = colors_set1[c(1:4)])
## Warning: Removed 441 rows containing missing values (geom_rect).

PlotDataTracksLight(input_tib = tib, 
                    exp_names = c("HCT116_AID_RO_Ki67",
                                  "HCT116_AID_RO_LMNB1",
                                  "HCT116_AID_RO_IAA_LMNB1",
                                  "HCT116_AID_RO_H3K27me3",
                                  "HCT116_AID_RO_IAA_H3K27me3",
                                  "HCT116_AID_RO_H3K9me3",
                                  "HCT116_AID_RO_IAA_H3K9me3"),
                    centromeres,
                    color_groups = c(HCT116_AID_RO_Ki67 = "Ki67",
                                     HCT116_AID_RO_LMNB1 = "LMNB1",
                                     HCT116_AID_RO_IAA_LMNB1 = "LMNB1",
                                     HCT116_AID_RO_H3K27me3 = "H3K27me3",
                                     HCT116_AID_RO_IAA_H3K27me3 = "H3K27me3",
                                     HCT116_AID_RO_H3K9me3 = "H3K9me3",
                                     HCT116_AID_RO_IAA_H3K9me3 = "H3K9me3"),
                    smooth = 9, plot_chr = "chr11", fix_yaxis = T,
                    plot_start = 25e6, plot_end = 100e6,
                    color_list = colors_set1[c(1:4)])
## Warning: Removed 395 rows containing missing values (geom_rect).

# Faceted plots
PlotDataTracksLightFaceted(input_tib = tib, 
                           exp_names = c("HCT116_AID_doublethy_Ki67",
                                         "HCT116_AID_doublethy_LMNB1",
                                         "HCT116_AID_doublethy_IAA_LMNB1",
                                         "HCT116_AID_doublethy_H3K27me3",
                                         "HCT116_AID_doublethy_IAA_H3K27me3",
                                         "HCT116_AID_doublethy_H3K9me3",
                                         "HCT116_AID_doublethy_IAA_H3K9me3"),
                           centromeres,
                           color_groups = c(HCT116_AID_doublethy_Ki67 = "minusIAA",
                                            HCT116_AID_doublethy_LMNB1 = "minusIAA",
                                            HCT116_AID_doublethy_IAA_LMNB1 = "plusIAA",
                                            HCT116_AID_doublethy_H3K27me3 = "minusIAA",
                                            HCT116_AID_doublethy_IAA_H3K27me3 = "plusIAA",
                                            HCT116_AID_doublethy_H3K9me3 = "minusIAA",
                                            HCT116_AID_doublethy_IAA_H3K9me3 = "plusIAA"),
                           smooth = 9, plot_chr = "chr3", fix_yaxis = F,
                           plot_start = 40e6, plot_end = 91e6,
                           color_list = c("grey20", "red"),
                           facet_groups = c(HCT116_AID_doublethy_Ki67 = "Ki67",
                                            HCT116_AID_doublethy_LMNB1 = "LMNB1",
                                            HCT116_AID_doublethy_IAA_LMNB1 = "LMNB1",
                                            HCT116_AID_doublethy_H3K27me3 = "H3K27me3",
                                            HCT116_AID_doublethy_IAA_H3K27me3 = "H3K27me3",
                                            HCT116_AID_doublethy_H3K9me3 = "H3K9me3",
                                            HCT116_AID_doublethy_IAA_H3K9me3 = "H3K9me3"))
## Warning: Removed 5 row(s) containing missing values (geom_path).

PlotDataTracksLightFaceted(input_tib = tib, 
                           exp_names = c("HCT116_AID_RO_Ki67",
                                         "HCT116_AID_RO_LMNB1",
                                         "HCT116_AID_RO_IAA_LMNB1",
                                         "HCT116_AID_RO_H3K27me3",
                                         "HCT116_AID_RO_IAA_H3K27me3",
                                         "HCT116_AID_RO_H3K9me3",
                                         "HCT116_AID_RO_IAA_H3K9me3"),
                           centromeres,
                           color_groups = c(HCT116_AID_RO_Ki67 = "minusIAA",
                                            HCT116_AID_RO_LMNB1 = "minusIAA",
                                            HCT116_AID_RO_IAA_LMNB1 = "plusIAA",
                                            HCT116_AID_RO_H3K27me3 = "minusIAA",
                                            HCT116_AID_RO_IAA_H3K27me3 = "plusIAA",
                                            HCT116_AID_RO_H3K9me3 = "minusIAA",
                                            HCT116_AID_RO_IAA_H3K9me3 = "plusIAA"),
                           smooth = 9, plot_chr = "chr3", fix_yaxis = F,
                           plot_start = 40e6, plot_end = 91e6,
                           color_list = c("grey20", "red"),
                           facet_groups = c(HCT116_AID_RO_Ki67 = "Ki67",
                                            HCT116_AID_RO_LMNB1 = "LMNB1",
                                            HCT116_AID_RO_IAA_LMNB1 = "LMNB1",
                                            HCT116_AID_RO_H3K27me3 = "H3K27me3",
                                            HCT116_AID_RO_IAA_H3K27me3 = "H3K27me3",
                                            HCT116_AID_RO_H3K9me3 = "H3K9me3",
                                            HCT116_AID_RO_IAA_H3K9me3 = "H3K9me3"))
## Warning: Removed 5 row(s) containing missing values (geom_path).

# Only the difference
PlotDataTracksLightFaceted(input_tib = tib, 
                           exp_names = c("doublethy_LMNB1_diff",
                                         "doublethy_H3K27me3_diff",
                                         "doublethy_H3K9me3_diff"),
                           centromeres,
                           color_groups = c(doublethy_LMNB1_diff = "LMNB1",
                                            doublethy_H3K27me3_diff = "H3K27me3",
                                            doublethy_H3K9me3_diff = "H3K9me3"),
                           smooth = 51, plot_chr = "chr3", 
                           plot_start = 40e6, plot_end = 91e6,
                           color_list = colors_set1[c(2:4)],
                           facet_groups = c(doublethy_LMNB1_diff = "diff",
                                            doublethy_H3K27me3_diff = "diff",
                                            doublethy_H3K9me3_diff = "diff"))

PlotDataTracksLightFaceted(input_tib = tib, 
                           exp_names = c("RO_LMNB1_diff",
                                         "RO_H3K27me3_diff",
                                         "RO_H3K9me3_diff"),
                           centromeres,
                           color_groups = c(RO_LMNB1_diff = "LMNB1",
                                            RO_H3K27me3_diff = "H3K27me3",
                                            RO_H3K9me3_diff = "H3K9me3"),
                           smooth = 51, plot_chr = "chr3", 
                           plot_start = 40e6, plot_end = 91e6,
                           color_list = colors_set1[c(2:4)],
                           facet_groups = c(RO_LMNB1_diff = "diff",
                                            RO_H3K27me3_diff = "diff",
                                            RO_H3K9me3_diff = "diff"))

# Other site
PlotDataTracksLightFaceted(input_tib = tib, 
                           exp_names = c("HCT116_AID_doublethy_Ki67",
                                         "HCT116_AID_doublethy_LMNB1",
                                         "HCT116_AID_doublethy_IAA_LMNB1",
                                         "HCT116_AID_doublethy_H3K27me3",
                                         "HCT116_AID_doublethy_IAA_H3K27me3",
                                         "HCT116_AID_doublethy_H3K9me3",
                                         "HCT116_AID_doublethy_IAA_H3K9me3"),
                           centromeres,
                           color_groups = c(HCT116_AID_doublethy_Ki67 = "minusIAA",
                                            HCT116_AID_doublethy_LMNB1 = "minusIAA",
                                            HCT116_AID_doublethy_IAA_LMNB1 = "plusIAA",
                                            HCT116_AID_doublethy_H3K27me3 = "minusIAA",
                                            HCT116_AID_doublethy_IAA_H3K27me3 = "plusIAA",
                                            HCT116_AID_doublethy_H3K9me3 = "minusIAA",
                                            HCT116_AID_doublethy_IAA_H3K9me3 = "plusIAA"),
                           smooth = 9, plot_chr = "chr15", fix_yaxis = F,
                           plot_start = 18e6, plot_end = 80e6,
                           color_list = c("grey20", "red"),
                           facet_groups = c(HCT116_AID_doublethy_Ki67 = "Ki67",
                                            HCT116_AID_doublethy_LMNB1 = "LMNB1",
                                            HCT116_AID_doublethy_IAA_LMNB1 = "LMNB1",
                                            HCT116_AID_doublethy_H3K27me3 = "H3K27me3",
                                            HCT116_AID_doublethy_IAA_H3K27me3 = "H3K27me3",
                                            HCT116_AID_doublethy_H3K9me3 = "H3K9me3",
                                            HCT116_AID_doublethy_IAA_H3K9me3 = "H3K9me3"))
## Warning: Removed 64 row(s) containing missing values (geom_path).

PlotDataTracksLightFaceted(input_tib = tib, 
                           exp_names = c("HCT116_AID_RO_Ki67",
                                         "HCT116_AID_RO_LMNB1",
                                         "HCT116_AID_RO_IAA_LMNB1",
                                         "HCT116_AID_RO_H3K27me3",
                                         "HCT116_AID_RO_IAA_H3K27me3",
                                         "HCT116_AID_RO_H3K9me3",
                                         "HCT116_AID_RO_IAA_H3K9me3"),
                           centromeres,
                           color_groups = c(HCT116_AID_RO_Ki67 = "minusIAA",
                                            HCT116_AID_RO_LMNB1 = "minusIAA",
                                            HCT116_AID_RO_IAA_LMNB1 = "plusIAA",
                                            HCT116_AID_RO_H3K27me3 = "minusIAA",
                                            HCT116_AID_RO_IAA_H3K27me3 = "plusIAA",
                                            HCT116_AID_RO_H3K9me3 = "minusIAA",
                                            HCT116_AID_RO_IAA_H3K9me3 = "plusIAA"),
                           smooth = 9, plot_chr = "chr15", fix_yaxis = F,
                           plot_start = 18e6, plot_end = 80e6,
                           color_list = c("grey20", "red"),
                           facet_groups = c(HCT116_AID_RO_Ki67 = "Ki67",
                                            HCT116_AID_RO_LMNB1 = "LMNB1",
                                            HCT116_AID_RO_IAA_LMNB1 = "LMNB1",
                                            HCT116_AID_RO_H3K27me3 = "H3K27me3",
                                            HCT116_AID_RO_IAA_H3K27me3 = "H3K27me3",
                                            HCT116_AID_RO_H3K9me3 = "H3K9me3",
                                            HCT116_AID_RO_IAA_H3K9me3 = "H3K9me3"))
## Warning: Removed 32 row(s) containing missing values (geom_path).

PlotDataTracksLightFaceted(input_tib = tib %>%
                             mutate(delta_lmnb1 = HCT116_AID_doublethy_IAA_LMNB1 - 
                                      HCT116_AID_doublethy_LMNB1), 
                           exp_names = c("HCT116_AID_doublethy_Ki67",
                                         "HCT116_AID_doublethy_LMNB1",
                                         "HCT116_AID_doublethy_IAA_LMNB1",
                                         "delta_lmnb1"),
                           centromeres,
                           color_groups = c(HCT116_AID_doublethy_Ki67 = "minusIAA",
                                            HCT116_AID_doublethy_LMNB1 = "minusIAA",
                                            HCT116_AID_doublethy_IAA_LMNB1 = "plusIAA", 
                                            delta_lmnb1 = "delta"),
                           smooth = 51, plot_chr = "chr3", fix_yaxis = F,
                           #plot_start = 40e6, plot_end = 91e6,
                           color_list = c("grey20", "red", "grey50"),
                           facet_groups = c(HCT116_AID_doublethy_Ki67 = "Ki67",
                                            HCT116_AID_doublethy_LMNB1 = "LMNB1",
                                            HCT116_AID_doublethy_IAA_LMNB1 = "LMNB1",
                                            delta_lmnb1 = "delta_LMNB1"))

PlotDataTracksLightFaceted(input_tib = tib %>%
                             mutate(delta_lmnb1 = HCT116_AID_RO_IAA_LMNB1 - 
                                      HCT116_AID_RO_LMNB1), 
                           exp_names = c("HCT116_AID_RO_Ki67",
                                         "HCT116_AID_RO_LMNB1",
                                         "HCT116_AID_RO_IAA_LMNB1",
                                         "delta_lmnb1"),
                           centromeres,
                           color_groups = c(HCT116_AID_RO_Ki67 = "minusIAA",
                                            HCT116_AID_RO_LMNB1 = "minusIAA",
                                            HCT116_AID_RO_IAA_LMNB1 = "plusIAA",
                                            delta_lmnb1 = "delta"),
                           smooth = 51, plot_chr = "chr3", fix_yaxis = F,
                           #plot_start = 40e6, plot_end = 91e6,
                           color_list = c("grey20", "red", "grey50"),
                           facet_groups = c(HCT116_AID_RO_Ki67 = "Ki67",
                                            HCT116_AID_RO_LMNB1 = "LMNB1",
                                            HCT116_AID_RO_IAA_LMNB1 = "LMNB1",
                                            delta_lmnb1 = "delta_LMNB1"))

# Original data
PlotDataTracksLight(input_tib = tib, 
                    exp_names = c("HCT116_AID_doublethy_Ki67",
                                  "HCT116_AID_doublethy_LMNB1",
                                  "HCT116_AID_doublethy_IAA_LMNB1"),
                    centromeres,
                    color_groups = c(HCT116_AID_doublethy_Ki67 = "Ki67",
                                     HCT116_AID_doublethy_LMNB1 = "LMNB1",
                                     HCT116_AID_doublethy_IAA_LMNB1 = "LMNB1"),
                    smooth = 9, plot_chr = "chr3", fix_yaxis = T,
                    #plot_start = 25e6, plot_end = 100e6,
                    color_list = colors_set1[c(1:2)])
## Warning: Removed 116 rows containing missing values (geom_rect).

PlotDataTracksLight(input_tib = tib, 
                    exp_names = c("HCT116_AID_RO_Ki67",
                                  "HCT116_AID_RO_LMNB1",
                                  "HCT116_AID_RO_IAA_LMNB1"),
                    centromeres,
                    color_groups = c(HCT116_AID_RO_Ki67 = "Ki67",
                                     HCT116_AID_RO_LMNB1 = "LMNB1",
                                     HCT116_AID_RO_IAA_LMNB1 = "LMNB1"),
                    smooth = 9, plot_chr = "chr3", fix_yaxis = T,
                    #plot_start = 25e6, plot_end = 100e6,
                    color_list = colors_set1[c(1:2)])
## Warning: Removed 116 rows containing missing values (geom_rect).

for (chr in chromosomes) {
  # doublethy
  PlotDataTracksLight(input_tib = tib, 
                      exp_names = c("HCT116_AID_doublethy_Ki67",
                                    "HCT116_AID_doublethy_LMNB1",
                                    "doublethy_LMNB1_diff",
                                    "doublethy_H3K27me3_diff",
                                    "doublethy_H3K9me3_diff"),
                      centromeres, 
                      color_groups = c(HCT116_AID_doublethy_Ki67 = "Ki67",
                                       HCT116_AID_doublethy_LMNB1 = "LMNB1",
                                       doublethy_LMNB1_diff = "LMNB1",
                                       doublethy_H3K27me3_diff = "H3K27me3",
                                       doublethy_H3K9me3_diff = "H3K9me3"),
                      smooth = 15, plot_chr = chr, fix_yaxis = F,
                      color_list = colors_set1[c(1:4)])
  
  # RO
  PlotDataTracksLight(input_tib = tib, 
                      exp_names = c("HCT116_AID_RO_Ki67",
                                    "HCT116_AID_RO_LMNB1",
                                    "RO_LMNB1_diff",
                                    "RO_H3K27me3_diff",
                                    "RO_H3K9me3_diff"),
                      centromeres,
                      color_groups = c(HCT116_AID_RO_Ki67 = "Ki67",
                                       HCT116_AID_RO_LMNB1 = "LMNB1",
                                       RO_LMNB1_diff = "LMNB1",
                                       RO_H3K27me3_diff = "H3K27me3",
                                       RO_H3K9me3_diff = "H3K9me3"),
                      smooth = 15, plot_chr = chr, fix_yaxis = F,
                      color_list = colors_set1[c(1:4)])
}
## Warning: Removed 1740 rows containing missing values (geom_rect).

## Warning: Removed 1739 rows containing missing values (geom_rect).

## Warning: Removed 165 rows containing missing values (geom_rect).

## Warning: Removed 165 rows containing missing values (geom_rect).

## Warning: Removed 145 rows containing missing values (geom_rect).

## Warning: Removed 145 rows containing missing values (geom_rect).

## Warning: Removed 145 rows containing missing values (geom_rect).

## Warning: Removed 145 rows containing missing values (geom_rect).

## Warning: Removed 184 rows containing missing values (geom_rect).

## Warning: Removed 170 rows containing missing values (geom_rect).

## Warning: Removed 105 rows containing missing values (geom_rect).

## Warning: Removed 105 rows containing missing values (geom_rect).

## Warning: Removed 205 rows containing missing values (geom_rect).

## Warning: Removed 205 rows containing missing values (geom_rect).

## Warning: Removed 109 rows containing missing values (geom_rect).

## Warning: Removed 105 rows containing missing values (geom_rect).

## Warning: Removed 1670 rows containing missing values (geom_rect).

## Warning: Removed 1675 rows containing missing values (geom_rect).

## Warning: Removed 145 rows containing missing values (geom_rect).

## Warning: Removed 145 rows containing missing values (geom_rect).

## Warning: Removed 292 rows containing missing values (geom_rect).

## Warning: Removed 255 rows containing missing values (geom_rect).

## Warning: Removed 185 rows containing missing values (geom_rect).

## Warning: Removed 184 rows containing missing values (geom_rect).

## Warning: Removed 1640 rows containing missing values (geom_rect).

## Warning: Removed 1634 rows containing missing values (geom_rect).

## Warning: Removed 1702 rows containing missing values (geom_rect).

## Warning: Removed 1679 rows containing missing values (geom_rect).

## Warning: Removed 1845 rows containing missing values (geom_rect).

## Warning: Removed 1830 rows containing missing values (geom_rect).

## Warning: Removed 930 rows containing missing values (geom_rect).

## Warning: Removed 915 rows containing missing values (geom_rect).

## Warning: Removed 320 rows containing missing values (geom_rect).

## Warning: Removed 312 rows containing missing values (geom_rect).

## Warning: Removed 479 rows containing missing values (geom_rect).

## Warning: Removed 479 rows containing missing values (geom_rect).

## Warning: Removed 200 rows containing missing values (geom_rect).

## Warning: Removed 189 rows containing missing values (geom_rect).

## Warning: Removed 145 rows containing missing values (geom_rect).

## Warning: Removed 143 rows containing missing values (geom_rect).

## Warning: Removed 608 rows containing missing values (geom_rect).

## Warning: Removed 610 rows containing missing values (geom_rect).

## Warning: Removed 1172 rows containing missing values (geom_rect).

## Warning: Removed 1177 rows containing missing values (geom_rect).

## Warning: Removed 250 rows containing missing values (geom_rect).

## Warning: Removed 245 rows containing missing values (geom_rect).

2. Differences between S-phase and G2 in HCT116 cells

Can I reproduce the cell cycle effects for Ki67 interactions?

# Plot chromosomal differences
tib %>%
  group_by(seqnames) %>%
  dplyr::summarise(Ki67 = mean(RO_vs_doublethy_Ki67, na.rm = T)) %>%
  ungroup() %>%
  mutate(seqnames = factor(seqnames, levels = chrom_sizes$seqnames)) %>%
  ggplot(aes(x = seqnames, y = Ki67)) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  geom_point() +
  ylab("Difference G2 versus S-phase") +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

# How about for all marks?
tib %>%
  group_by(seqnames) %>%
  dplyr::summarise(Ki67 = mean(RO_vs_doublethy_Ki67, na.rm = T),
                   LMNB1 = mean(RO_vs_doublethy_LMNB1, na.rm = T),
                   H3K27me3 = mean(RO_vs_doublethy_H3K27me3, na.rm = T),
                   H3K9me3 = mean(RO_vs_doublethy_H3K9me3, na.rm = T)) %>%
  ungroup() %>%
  gather(target, value, -seqnames) %>%
  mutate(seqnames = factor(seqnames, chrom_sizes$seqnames),
         target = factor(target, levels(padamid_metadata_combined$target))) %>%
  ggplot(aes(x = seqnames, y = value, col = target)) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  geom_point() +
  ylab("Difference G2 versus S-phase") +
  scale_color_manual(values = colors_set1) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

Yes. That’s good.

3. Ki67 score versus auxin-induced differences

So, can I show that auxin depletion results in perturbed lamina interactions in regions that are high in initial Ki67 binding?

Some issues: * Data is quite noisy, in particular Ki67 tracks * Small chromosomes are high in Ki67, but seem to stay interior after Ki67 depletion as well

# Plot Ki67 score versus the differences?
PlotScatterBinned(tib, n1 = "HCT116_AID_doublethy_Ki67", n2 = "doublethy_LMNB1_diff", 
                  smooth = 1, count_range = c(0, 1000))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, n1 = "HCT116_AID_doublethy_Ki67", n2 = "doublethy_H3K27me3_diff", 
                  smooth = 1, count_range = c(0, 1000))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, n1 = "HCT116_AID_doublethy_Ki67", n2 = "doublethy_H3K9me3_diff", 
                  smooth = 1, count_range = c(0, 1000))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, n1 = "HCT116_AID_RO_Ki67", n2 = "RO_LMNB1_diff", 
                  smooth = 1, count_range = c(0, 1000))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, n1 = "HCT116_AID_RO_Ki67", n2 = "RO_H3K27me3_diff", 
                  smooth = 1, count_range = c(0, 1000))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, n1 = "HCT116_AID_RO_Ki67", n2 = "RO_H3K9me3_diff", 
                  smooth = 1, count_range = c(0, 1000))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# This seems to noisy

# Plot LaminB1 scores, highlight initial Ki67
PlotScatterBinned(tib, smooth = 1,
                  n1 = "HCT116_AID_doublethy_LMNB1", 
                  n2 = "HCT116_AID_doublethy_IAA_LMNB1",
                  color_by = tib$HCT116_AID_doublethy_Ki67, ylimits_col = c(-2.5, 2.5))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 1,
                  n1 = "HCT116_AID_doublethy_H3K27me3", 
                  n2 = "HCT116_AID_doublethy_IAA_H3K27me3",
                  color_by = tib$HCT116_AID_doublethy_Ki67, ylimits_col = c(-2.5, 2.5))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 1,
                  n1 = "HCT116_AID_doublethy_H3K9me3", 
                  n2 = "HCT116_AID_doublethy_IAA_H3K9me3",
                  color_by = tib$HCT116_AID_doublethy_Ki67, ylimits_col = c(-2.5, 2.5))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 1,
                  n1 = "HCT116_AID_RO_LMNB1", 
                  n2 = "HCT116_AID_RO_IAA_LMNB1",
                  color_by = tib$HCT116_AID_RO_Ki67, ylimits_col = c(-2.5, 2.5))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 1,
                  n1 = "HCT116_AID_RO_H3K27me3", 
                  n2 = "HCT116_AID_RO_IAA_H3K27me3",
                  color_by = tib$HCT116_AID_RO_Ki67, ylimits_col = c(-2.5, 2.5))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 1,
                  n1 = "HCT116_AID_RO_H3K9me3", 
                  n2 = "HCT116_AID_RO_IAA_H3K9me3",
                  color_by = tib$HCT116_AID_RO_Ki67, ylimits_col = c(-2.5, 2.5))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# Plot LaminB1 scores, highlight initial Ki67 - 3 bins smoothing?
PlotScatterBinned(tib, smooth = 3,
                  n1 = "HCT116_AID_doublethy_LMNB1", 
                  n2 = "HCT116_AID_doublethy_IAA_LMNB1",
                  color_by = tib$HCT116_AID_doublethy_Ki67, ylimits_col = c(-2.8, 2.8))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "HCT116_AID_doublethy_H3K27me3", 
                  n2 = "HCT116_AID_doublethy_IAA_H3K27me3",
                  color_by = tib$HCT116_AID_doublethy_Ki67, ylimits_col = c(-2.8, 2.8))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "HCT116_AID_doublethy_H3K9me3", 
                  n2 = "HCT116_AID_doublethy_IAA_H3K9me3",
                  color_by = tib$HCT116_AID_doublethy_Ki67, ylimits_col = c(-2.8, 2.8))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "HCT116_AID_RO_LMNB1", 
                  n2 = "HCT116_AID_RO_IAA_LMNB1",
                  color_by = tib$HCT116_AID_RO_Ki67, ylimits_col = c(-2.8, 2.8))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "HCT116_AID_RO_H3K27me3", 
                  n2 = "HCT116_AID_RO_IAA_H3K27me3",
                  color_by = tib$HCT116_AID_RO_Ki67, ylimits_col = c(-2.8, 2.8))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib, smooth = 3,
                  n1 = "HCT116_AID_RO_H3K9me3", 
                  n2 = "HCT116_AID_RO_IAA_H3K9me3",
                  color_by = tib$HCT116_AID_RO_Ki67, ylimits_col = c(-2.8, 2.8))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# Looks better

# Something else: correlation per chromosome with the difference
tib_cor <- tib %>%
  group_by(seqnames) %>%
  drop_na() %>%
  dplyr::summarise(doublethy_LMNB1 = cor(HCT116_AID_doublethy_Ki67,
                                         doublethy_LMNB1_diff),
                   doublethy_H3K27me3 = cor(HCT116_AID_doublethy_Ki67,
                                            doublethy_H3K27me3_diff),
                   doublethy_H3K9me3 = cor(HCT116_AID_doublethy_Ki67,
                                           doublethy_H3K9me3_diff),
                   RO_LMNB1 = cor(HCT116_AID_RO_Ki67,
                                  RO_LMNB1_diff),
                   RO_H3K27me3 = cor(HCT116_AID_RO_Ki67,
                                     RO_H3K27me3_diff),
                   RO_H3K9me3 = cor(HCT116_AID_RO_Ki67,
                                    RO_H3K9me3_diff)) %>%
  ungroup() %>%
  gather(key, value, -seqnames) %>%
  separate(key, c("condition", "target"), remove = F) %>%
  mutate(condition = factor(condition, levels(padamid_metadata_combined$condition)),
         target = factor(target, levels(padamid_metadata_combined$target)),
         seqnames = factor(seqnames, chrom_sizes$seqnames))

tib_cor %>%
  ggplot(aes(x = seqnames, y = value, col = target)) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  geom_point(size = 3) +
  facet_grid(. ~ condition) +
  ylab("Pearson correlation between Ki67 and auxin-induced difference") +
  scale_color_manual(values = colors_set1[2:4]) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

# Mweh. Difficult to interpret.

# How about boxplots per Ki67 deciles?
tib_decile_doublethy <- tib %>%
    drop_na() %>%
    mutate(ki67_decile = case_when(HCT116_AID_doublethy_Ki67 > 1 ~ "high",
                                   HCT116_AID_doublethy_Ki67 < -1 ~ "low",
                                   T ~ "medium"),
           ki67_decile = factor(ki67_decile, c("low", "medium", "high"))) %>%
    dplyr::select(seqnames, ki67_decile, starts_with("doublethy")) %>%
    gather(key, value, -seqnames, -ki67_decile) %>%
    mutate(key = str_remove(key, "_diff"),
           key = str_remove(key, ".*_"),
           key = factor(key, levels(padamid_metadata_combined$target)))

tib_decile_doublethy %>%
    ggplot(aes(x = ki67_decile, y = value, fill = ki67_decile)) +
    geom_boxplot(outlier.shape = NA) +
    coord_cartesian(ylim = c(-1.3, 1.3)) +
    facet_grid(. ~ key) +
    theme_bw() +
    theme(aspect.ratio = 1)

tib_decile_RO <- tib %>%
    drop_na() %>%
    arrange(HCT116_AID_RO_Ki67) %>%
    mutate(ki67_decile = case_when(HCT116_AID_RO_Ki67 > 1 ~ "high",
                                   HCT116_AID_RO_Ki67 < -1 ~ "low",
                                   T ~ "medium"),
           ki67_decile = factor(ki67_decile, c("low", "medium", "high"))) %>%
    dplyr::select(seqnames, ki67_decile, starts_with("RO") & ends_with("diff")) %>%
    gather(key, value, -seqnames, -ki67_decile) %>%
    mutate(key = str_remove(key, "_diff"),
           key = str_remove(key, ".*_"),
           key = factor(key, levels(padamid_metadata_combined$target)))

tib_decile_RO %>%
    ggplot(aes(x = ki67_decile, y = value, fill = ki67_decile)) +
    geom_boxplot(outlier.shape = NA) +
    coord_cartesian(ylim = c(-1.3, 1.3)) +
    facet_grid(. ~ key) +
    theme_bw() +
    theme(aspect.ratio = 1)

# Mweh. Not clear.

# Something different even - correlations not by chromosome
tib_cor <- tib %>%
  # Smoothing of 3 bins?
  mutate_at(vars(contains("vs"), contains("diff")), 
              function(x) runmean(x, k = 3)) %>%
  drop_na() %>%
  dplyr::summarise(doublethy_LMNB1 = cor(HCT116_AID_doublethy_Ki67,
                                         doublethy_LMNB1_diff),
                   doublethy_H3K27me3 = cor(HCT116_AID_doublethy_Ki67,
                                            doublethy_H3K27me3_diff),
                   doublethy_H3K9me3 = cor(HCT116_AID_doublethy_Ki67,
                                           doublethy_H3K9me3_diff),
                   RO_LMNB1 = cor(HCT116_AID_RO_Ki67,
                                  RO_LMNB1_diff),
                   RO_H3K27me3 = cor(HCT116_AID_RO_Ki67,
                                     RO_H3K27me3_diff),
                   RO_H3K9me3 = cor(HCT116_AID_RO_Ki67,
                                    RO_H3K9me3_diff)) %>%
  gather(key, value) %>%
  separate(key, c("condition", "target"), remove = F) %>%
  mutate(condition = factor(condition, levels(padamid_metadata_combined$condition)),
         target = factor(target, levels(padamid_metadata_combined$target)))

tib_cor %>%
  ggplot(aes(x = condition, y = value, col = target)) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  geom_point(size = 3) +
  xlab("") +
  ylab("Pearson correlation between Ki67 and auxin-induced difference") +
  scale_color_manual(values = colors_set1[2:4]) +
  theme_bw() +
  theme(aspect.ratio = 1.5,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

tib_cor %>%
  print(n = 20)
## # A tibble: 6 Ă— 4
##   key                condition target     value
##   <chr>              <fct>     <fct>      <dbl>
## 1 doublethy_LMNB1    doublethy LMNB1     0.291 
## 2 doublethy_H3K27me3 doublethy H3K27me3  0.246 
## 3 doublethy_H3K9me3  doublethy H3K9me3  -0.0349
## 4 RO_LMNB1           RO        LMNB1     0.205 
## 5 RO_H3K27me3        RO        H3K27me3  0.127 
## 6 RO_H3K9me3         RO        H3K9me3   0.124
# Easy to explain, easy to interpret, show only a fraction of the data and the
# entire story
#
# Let's go for it




# Include distance to centromere / smoothing of Ki-67 differences
tib_tmp <- tib %>%
  mutate(distance_to_centromere = mcols(distanceToNearest(as(tib, "GRanges"),
                                                          centromeres))$distance,
         distance_to_centromere = distance_to_centromere / 1e6) %>%
  mutate_at(vars(matches("HCT116|diff")), function(x) runmean(x, k = 9))

# Scatterplot
PlotScatterBinned(tib_tmp, smooth = 3,
                  n1 = "HCT116_AID_doublethy_Ki67",
                  n2 = "doublethy_LMNB1_diff", 
                  color_by = log10(tib_tmp$distance_to_centromere+1), ylimits_col = c(-5, 5))
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

# All scatterplots
tib_cor <- tibble()
conditions <- c("doublethy", "RO")
targets <- c("LMNB1", "H3K27me3", "H3K9me3")

for (condition in conditions) {
  for (target in targets) {
    
    # Get the data
    tib_tmp <- tib_padamid_combined %>%
      dplyr::select(1:3,
                    matches("AID") &
                      matches(condition) &
                      matches(target)) %>%
      rename_at(4:5, ~ c("DMSO", "IAA")) %>%
      drop_na()
    
    # Scatterplot
    #PlotScatterBinned(tib_tmp, "DMSO", "IAA", count_range = c(0, 750))
    
    # Determine pearson correlation
    cor_tmp <- cor(tib_tmp$DMSO, tib_tmp$IAA)
    
    # Add this to tib_cor object
    tib_cor <- bind_rows(tib_cor,
                         tibble(condition = condition,
                                target = target, 
                                cor = cor_tmp))
  }
}

# Add factors
tib_cor <- tib_cor %>%
  mutate(condition = factor(condition, levels = conditions),
         target = factor(target, levels = targets))

# Combined plot
tib_cor %>%
  ggplot(aes(x = target, y = cor, col = condition)) +
  geom_point(size = 3) +
  xlab("") +
  ylab("Pearson correlation") +
  scale_color_manual(values = colors_set3) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 1)) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

4. Near centromeric summaries

# Add distance to centromeres
dis <- distanceToNearest(as(tib, "GRanges"), centromeres)
tib$distance_to_centromere <- mcols(dis)$distance / 1e6

# Summary per chromosome
tib_summary <- tib %>%
  filter(distance_to_centromere < 2) %>%
  drop_na() %>%
  gather(key, value, contains("diff")) %>%
  group_by(seqnames, key) %>%
  summarise(mean = mean(value)) %>%
  ungroup() %>%
  separate(key, c("condition", "target"), remove = F) %>%
  mutate(condition = factor(condition, c("doublethy", "RO")),
         target = factor(target, levels(padamid_metadata_combined$target)),
         seqnames = factor(seqnames, chrom_sizes$seqnames),
         rdna = seqnames %in% paste0("chr", c(13, 14, 15, 21, 22))) %>%
  mutate(size = chrom_sizes$length[match(seqnames, 
                                         chrom_sizes$seqnames)],
         size = size / 1e6)
## `summarise()` has grouped output by 'seqnames'. You can override using the `.groups` argument.
## Warning: Expected 2 pieces. Additional pieces discarded in 138 rows [1, 2, 3, 4,
## 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].
# Plot
tib_summary %>%
  ggplot(aes(x = seqnames, y = mean, col = target, shape = rdna)) +
  geom_point(size = 2) +
  geom_hline(yintercept = 0, col = "black") +
  facet_grid(. ~ condition) +
  xlab("") +
  ylab("Ki67 score near centromeres") +
  scale_color_manual(values = colors_set1[2:4]) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

Is this due to DNA accessibility changes?

# Load Dam counts, for some samples
# List files
count_files <- dir("results/tracks/counts/bin-50kb/", pattern = "Dam", 
                   recursive = T, full.names = T)
count_files <- grep("AID", count_files, value = T)
count_files <- grep("doublethy|RO", count_files, value = T)


# Prepare into metadata
count_metadata <- tibble(file = count_files) %>%
  mutate(sample = basename(str_remove(file, "\\..*")),
         sample = str_remove(sample, "_Dam.*"),
         sample = str_remove(sample, ".*AID_")) %>%
  mutate(condition = case_when(str_detect(sample, "double") ~ "doublethy",
                               T ~ "RO"),
         treatment = case_when(str_detect(sample, "IAA") ~ "plusIAA",
                               T ~ "minusIAA"),
         replicate = str_remove(sample, ".*_"),
         class = paste(condition, treatment, sep = "_"))

# Load count bigwig files - combine into one tibble
tmp <- tibble(seqnames = character(),
              start = numeric(),
              end = numeric())

for (i in 1:nrow(count_metadata)) {
  f_name <- count_metadata$sample[i]
  f_tib <- as_tibble(import(count_metadata$file[i])) %>%
    dplyr::select(-width, -strand) %>%
    rename(score = f_name) %>%
    mutate(seqnames = as.character(seqnames))
  tmp <- full_join(tmp, f_tib)
}
## Joining, by = c("seqnames", "start", "end")
## Joining, by = c("seqnames", "start", "end")
## Joining, by = c("seqnames", "start", "end")
## Joining, by = c("seqnames", "start", "end")
## Joining, by = c("seqnames", "start", "end")
## Joining, by = c("seqnames", "start", "end")
## Joining, by = c("seqnames", "start", "end")
## Joining, by = c("seqnames", "start", "end")
## Joining, by = c("seqnames", "start", "end")
## Joining, by = c("seqnames", "start", "end")
## Joining, by = c("seqnames", "start", "end")
## Joining, by = c("seqnames", "start", "end")
# Mean of replicates and bigger bins
ovl <- findOverlaps(as(tmp, "GRanges"), as(tib, "GRanges"), 
                    select = "arbitrary")
tmp <- tmp %>%
  add_column(bin = ovl) %>%
  gather(key, value, count_metadata$sample) %>%
  mutate(class = count_metadata$class[match(key, count_metadata$sample)]) %>%
  group_by(bin, class) %>%
  dplyr::summarise(mean = mean(value, na.rm = T)) %>%
  spread(class, mean) %>%
  arrange(bin) %>%
  drop_na(bin)
## `summarise()` has grouped output by 'bin'. You can override using the `.groups` argument.
# Prepare count GRanges
gr_count <- as(tib, "GRanges")
mcols(gr_count) <- NULL
mcols(gr_count)[, names(tmp)[2:ncol(tmp)]] <- NA
mcols(gr_count)[tmp$bin, names(tmp)[2:ncol(tmp)]] <- tmp[, 2:ncol(tmp)]

# Prepare count tibble & determine differences
tib_count <- as_tibble(gr_count) %>%
  dplyr::select(-strand, -width) %>%
  mutate(doublethy_diff = log2((doublethy_plusIAA+1) / (doublethy_minusIAA+1)),
         RO_diff = log2((RO_plusIAA+1) / (RO_minusIAA+1)))





# Repeat centromere exercise
# Add distance to centromeres
dis <- distanceToNearest(as(tib_count, "GRanges"), centromeres)
tib_count$distance_to_centromere <- mcols(dis)$distance / 1e6

# Summary per chromosome
tib_summary <- tib_count %>%
  filter(distance_to_centromere < 2) %>%
  drop_na() %>%
  gather(key, value, contains("diff")) %>%
  group_by(seqnames, key) %>%
  summarise(mean = mean(value)) %>%
  ungroup() %>%
  separate(key, c("condition"), remove = F) %>%
  mutate(condition = factor(condition, c("doublethy", "RO")),
         seqnames = factor(seqnames, chrom_sizes$seqnames),
         rdna = seqnames %in% paste0("chr", c(13, 14, 15, 21, 22))) %>%
  mutate(size = chrom_sizes$length[match(seqnames, 
                                         chrom_sizes$seqnames)],
         size = size / 1e6)
## `summarise()` has grouped output by 'seqnames'. You can override using the `.groups` argument.
## Warning: Expected 1 pieces. Additional pieces discarded in 46 rows [1, 2, 3, 4,
## 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].
# Plot
tib_summary %>%
  ggplot(aes(x = seqnames, y = mean, col = condition, shape = rdna)) +
  geom_point(size = 2) +
  geom_hline(yintercept = 0, col = "black") +
  xlab("") +
  ylab("Dam only differences near centromere") +
  scale_color_manual(values = colors_set3) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

# Alternative, correlate Ki67 versus the change in Dam signal
tib_combined <- full_join(tib, tib_count)
## Joining, by = c("seqnames", "start", "end", "distance_to_centromere")
PlotScatterBinned(tib_combined, "HCT116_AID_doublethy_Ki67", "doublethy_diff")
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib_combined, "HCT116_AID_RO_Ki67", "RO_diff")
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib_combined, "doublethy_minusIAA", "doublethy_plusIAA",
                  color_by = tib_combined$HCT116_AID_doublethy_Ki67,
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

PlotScatterBinned(tib_combined, "RO_minusIAA", "RO_plusIAA",
                  color_by = tib_combined$HCT116_AID_RO_Ki67,
                  remove_outliers = T)
## `summarise()` has grouped output by 'n1_bin'. You can override using the `.groups` argument.

Conclusions

Yes, Ki67 depletion results in increased LaminB1 interactions for regions initially high in Ki67. The differences are quite small, which means that this analysis is strongly affected by noisiness. I decided not to smooth anything.

I believe that the observations that I will present are legit, but might only tell a fraction of the complete story. Let’s go for it, though.

Session info

sessionInfo()
## R version 4.0.5 (2021-03-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] knitr_1.33           caTools_1.18.2       corrr_0.4.3         
##  [4] GGally_2.1.2         RColorBrewer_1.1-2   rtracklayer_1.50.0  
##  [7] GenomicRanges_1.42.0 GenomeInfoDb_1.26.7  IRanges_2.24.1      
## [10] S4Vectors_0.28.1     BiocGenerics_0.36.1  forcats_0.5.1       
## [13] stringr_1.4.0        dplyr_1.0.7          purrr_0.3.4         
## [16] readr_2.0.0          tidyr_1.1.3          tibble_3.1.6        
## [19] ggplot2_3.3.5        tidyverse_1.3.1     
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-7                matrixStats_0.60.0         
##  [3] fs_1.5.0                    lubridate_1.7.10           
##  [5] httr_1.4.2                  tools_4.0.5                
##  [7] backports_1.2.1             bslib_0.2.5.1              
##  [9] utf8_1.2.2                  R6_2.5.1                   
## [11] DBI_1.1.1                   colorspace_2.0-2           
## [13] withr_2.4.2                 tidyselect_1.1.1           
## [15] compiler_4.0.5              cli_3.1.0                  
## [17] rvest_1.0.1                 Biobase_2.50.0             
## [19] xml2_1.3.2                  DelayedArray_0.16.3        
## [21] labeling_0.4.2              sass_0.4.0                 
## [23] scales_1.1.1                digest_0.6.28              
## [25] Rsamtools_2.6.0             rmarkdown_2.11             
## [27] XVector_0.30.0              pkgconfig_2.0.3            
## [29] htmltools_0.5.1.1           MatrixGenerics_1.2.1       
## [31] highr_0.9                   dbplyr_2.1.1               
## [33] rlang_0.4.12                readxl_1.3.1               
## [35] rstudioapi_0.13             farver_2.1.0               
## [37] jquerylib_0.1.4             generics_0.1.0             
## [39] jsonlite_1.7.2              BiocParallel_1.24.1        
## [41] RCurl_1.98-1.3              magrittr_2.0.1             
## [43] GenomeInfoDbData_1.2.4      Matrix_1.3-4               
## [45] Rcpp_1.0.7                  munsell_0.5.0              
## [47] fansi_0.5.0                 lifecycle_1.0.1            
## [49] stringi_1.7.3               yaml_2.2.1                 
## [51] SummarizedExperiment_1.20.0 zlibbioc_1.36.0            
## [53] plyr_1.8.6                  grid_4.0.5                 
## [55] crayon_1.4.2                lattice_0.20-44            
## [57] Biostrings_2.58.0           haven_2.4.1                
## [59] hms_1.1.0                   pillar_1.6.4               
## [61] codetools_0.2-18            reprex_2.0.0               
## [63] XML_3.99-0.6                glue_1.5.0                 
## [65] evaluate_0.14               modelr_0.1.8               
## [67] vctrs_0.3.8                 tzdb_0.1.2                 
## [69] cellranger_1.1.0            gtable_0.3.0               
## [71] reshape_0.8.8               assertthat_0.2.1           
## [73] xfun_0.24                   broom_0.7.9                
## [75] GenomicAlignments_1.26.0    ellipsis_0.3.2